%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

% SET THE BELOW LINE TO THE OSL DIRECTORY:
%tilde='/home/mwoolrich';
tilde='/Users/woolrich';

osldir=[tilde '/homedir/matlab/osl1.2.beta.15'];
addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
% This specifies where our data is stored, what the data filenames are, and
% a few parameters for preprocessing such as how we will epoch the data.

datadir=[tilde '/homedir/matlab/osl_testdata_dir/faces_group_data_new_full']; % this is the directory the oat files will be stored in
spmfilesdir=[datadir '/spmfiles']; % this is the directory the SPM files will be stored in
structuralsdir=[datadir '/structurals']; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' datadir]; unix(cmd); % make dir to put the results in

clear fif_files spm_files structural_files

% Alternatively you can leave fif_files empty and set up a list of SPM MEEG
% object files:
is=[15:28];
spm_files=[];
for i=1:length(is),    
    spm_files{i}=[spmfilesdir '/espm8_meg' num2str(is(i)) '.mat'];
end;

% structural files:
sfs={'F9'    'F11'    'F15'    'F16'    'F20'    'F22'    'F25'    'M3'    'M6'    'M7'    'M12'    'M14'    'M18'    'M19'    'F4'    'F7'    'F10'    'F13'    'F18'    'F24'    'F26'    'M2'    'M4'    'M8'    'M10'    'M13'   'M17'    'M21'    'F3'    'F6'    'F12'    'F14'    'F17'    'F19'    'F23'    'M5'    'M9'    'M11'    'M15'    'M20'    'M22'    'M23'};
sfs{21}=''; % structural is missing
sfs{14}=''; % structural is missing
sfs{6}=''; % structural is missing
for i=1:length(is),
    if ~strcmp(sfs{is(i)},'')
        structural_files{i}=[structuralsdir '/' sfs{is(i)} '/struct.nii'];
    else
        structural_files{i}='';
    end;
end;

%%%%%%%%%%%%%%%%%%%
%% DO REGISTRATION AND RUN FORWARD MODEL BASED ON STRUCTURAL SCANS
% Before running the beamformer we need to compute the forward model for
% each subject based on their structural scan.
% Make sure you check the results look reasonable!

for i=1:length(spm_files),
    
        S2=[];
                
        S2.fid_label.nasion='Nasion';        
        S2.fid_label.lpa='LPA';
        S2.fid_label.rpa='RPA';
          
        S2.D = spm_files{i};    % requires .mat extension    
        S2.mri=structural_files{i}; % set S2.sMRI=''; if there is no structural available        
        S2.useheadshape=1;
                
        D=osl_forward_model(S2);
        
        D=osl_neuromag_grad_baseline_correction(spm_files{i});

end;
        
ca;
D=spm_eeg_load(spm_files{1});
%spm_eeg_inv_checkmeshes(D);
spm_eeg_inv_checkdatareg(D);
%spm_eeg_inv_checkforward(D, 1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP BEAMFORMER AND FIRST-LEVEL GLM OAT
% In this section we will do a wholebrain beamformer, followed by a trial-wise
% GLM that will correspond to a comparison of the ERFs for the different
% conditions.

%% SETUP THE OAT:
oat=[];
oat.source_recon.D_continuous=[]; % do not have continuous files
oat.source_recon.D_epoched=spm_files; % only epoched files
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.freq_range=[1 40]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.3];
oat.source_recon.method='beamform';
oat.source_recon.gridstep=8; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[datadir '/faces_group_norm'];
oat.source_recon.forward_meg='Single Sphere';
oat.source_recon.sessions_to_do=1:3;

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
%oat.first_level.contrast{4}=[0 1 0 0]'; % neutral
%oat.first_level.contrast{5}=[0 0 0 1]'; % fearful
%oat.first_level.contrast{6}=[0 -1 0 1]'; % fearful-neutral
oat.first_level.contrast_name={};
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
%oat.first_level.contrast_name{4}='neutral';
%oat.first_level.contrast_name{5}='fearful';
%oat.first_level.contrast_name{6}='fearful-neutral';

oat = osl_check_oat(oat);

%% run OAT
oat.source_recon.subjects_to_do=1:length(spm_files);
oat.first_level.subjects_to_do=oat.source_recon.sessions_to_do;

% for i=1:14, oat.first_level.results_fnames{i}=['subject' num2str(i) '_first_level']; end;
oat.to_do=[1 1 1 0];

oat = osl_run_oat(oat);

if(0),
    %%%%%%%%%%%%%%%%%%%
    %% OUTPUT SUBJECT'S NIFTII FILES
    % Having run the GLM on our source space data, we would like to inspect the
    % results for our single subject. 
    % We can do this by saving the contrast of parameter estimates (COPEs) and 
    % t-statistics for each of our contrasts to NIFTI images.

    S2=[];
    S2.oat=oat;
    S2.stats_fname=oat.first_level.results_fnames{1};
    S2.first_level_contrasts=[1,3]; % list of first level contrasts to output

    [statsdir times]=osl_save_nii_stats(S2);
end;

%%%%%%%%%%%%%%%%%%%
%% OUTPUT GROUP'S NIFTII FILES

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.first_level_contrasts=[3]; % list of first level contrasts to output

[statsdir,times]=osl_save_nii_stats(S2);
